rng(123)
global M theta delta Me Mc sigma X B k kappa r
Ind = 1000;          %Number of individuals in a city
T = 300;                %periods

M=zeros(Ind,T);        %Initial cash usage
B=zeros(Ind,T);        %benifits of switching
x=zeros(Ind,T);        %matrix of cash choice
Mc=1;
%lamda=0.8;   %network parameter



%%Parameters
r=(-log(0.9))/300;                 
theta=(-log(1-0.9))/70;
sigma=0.16;
delta=1/300;


% generate age distribution in the city
age1=randi(15,[round(Ind*old_ratio),1])+65;
age2=randi(44,[round(Ind*(1-old_ratio)),1])+20;
age=[age1;age2];




M(:,1)=Mc*ones(Ind,1)+normrnd(0,1,Ind,1); 
B(:,1)=(delta/(1-exp(-r*delta)))*(Me-Mc)*ones(Ind,1)++lamda*(1-old_ratio)+(((delta*exp(-theta*delta)))/(1-exp(-(r+theta)*delta)))*(Mc*ones(Ind,1)-M(:,1));
kappa_initial=prctile(B(:,1),old_ratio*100);
kappa=kappa_initial*ones(Ind,1);


for i=1:1000
    if B(i,1)>=kappa_initial
        x(i,1)=0;
    else
        x(i,1)=1;
    end
end

a=(1-exp(-theta*delta))*Mc;



for iteration=1:100
    iteration
for t = 2: T
    if t>=T/2 && t<=T/2+70
        shock=1;
    else
        shock=0;
    end
        

    
        M0=M(:,t-1);  
        B0=B(:,t-1);
     
        xm1=mean(x(:,t-1));
        R=normrnd(0,1,Ind,1);

      
   
      
      for j = 1: Ind
      
      M(j,t)=(1-exp(-theta*delta))*Mc+exp(-(theta)*delta)*M(j,t-1)+(delta^1/2)*sigma*R(j,1)-shock*shock_val*Me*old_ratio; 
      B(j,t)=(delta/(1-exp(-r*delta)))*(Me-Mc)+lamda*(1-xm1)+(((delta*exp(-theta*delta)))/(1-exp(-(r+theta)*delta)))*(Mc-M(j,t));
        
         if B(j,t)>=kappa(j,1) && x(j,t-1)>0
            x(j,t)=0;
         elseif B(j,t)>=kappa(j,1) && x(j,t-1)==0
            x(j,t)=0;
         elseif B(j,t)<kappa(j,1) && x(j,t-1)>0
            x(j,t)=1;
         else
            x(j,t)=0;
          end
      end
         
          eval(['B_mcc',num2str(iteration),'=B;']);
          eval(['M_mcc',num2str(iteration),'=M;']);
          eval(['x_mcc',num2str(iteration),'=x;']);
        
end
end

B=B_mcc1;
M=M_mcc1;
x=x_mcc1;
for i=2:iteration
eval(['B=B+B_mcc',num2str(i),';']);
eval(['M=M+M_mcc',num2str(i),';']);
eval(['x=x+x_mcc',num2str(i),';']);
end

B=B/iteration;
M=M/iteration;
x=x/iteration;

set(gcf,'color','w')
subplot(2,2,1);
M_mean=mean(M,1);
plot(M_mean(1,:),'-bo','Linewidth',1,'MarkerSize',1)
title('\sl{M_t}')   

subplot(2,2,2);
B_mean=mean(B,1);
plot(B_mean(1,:),'-bo','Linewidth',1,'MarkerSize',1)
title('\sl{B_t}')   

subplot(2,2,3);
kappa_mean=ones(1,300)*kappa_initial;
plot(kappa_mean(1,:),'-bo','Linewidth',1,'MarkerSize',1)
title('\kappa')   

subplot(2,2,4);
x_mean=mean(x,1);
plot(x_mean(1,:),'-bo','Linewidth',1,'MarkerSize',1)
title('\sl{x_t}')   
      